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SECTION  I 


INTRODUCTION 

AiResearch,  under  contract  to  USAF/APL  (Contract  Number 
F33615-76-C-2071 ) , has  been  developing  the  technology  to  perform 
transonic  analysis  of  compressors  with  and  without  splitter  vanes. 
This  has  generally  followed  five  phases  of  development,  proceed- 
ing from  a single-bladed  cascade  to  splittered  cascades,  and 
finally  to  a fully-rotating  compressor.  The  following  sections 
describe  the  basic  method  (description  of  geometry,  radial 
equilibrium,  and  relaxation  portions)  and  present  a comparison 
to  data. 

SECTION  II 
THE  METHOD 

By  necessity,  the  design  of  transonic  compressors  has  al- 
ways been  performed  utilizing  methods  that  are  incapable  of 
predicting  the  real  flow  fields  internal  to  such  devices. 
Consequently , there  has  been  a heavy  reliance  on  test-based 
development  programs.  The  cost  of  testing  has,  however, 
increased  so  rapidly  that  the  amount  of  testing,  and  therefore 
the  performance  of  new  components,  is  largely  controlled  by 
fiscal  constraints.  The  bright  spot,  in  what  would  otherwise 
be  a gloomy  outlook,  is  that  the  cost  of  computing  has  steadily 
declined  as  larger  and  larger  computers  are  developed.  This, 
coupled  with  vastly  improved  numerical  methods,  has  resulted, 
for  at  least  some  turbomachinery  devices,  in  the  replacement  of 
testing  by  computational  experiments.  However,  of  all  turbo- 
machinery devices,  the  transonic  compressor  is  one  of  the  most 
complicated.  Computational  developments  for  this  device 
initially  concentrated  on  transoni^  flow.  Among  these  develop- 
ments has  been  the  work  of  Dodge1 ' “ applying  transonic  relax- 
ation to  internal  flows. 

Rae^/^  has  applied  similar  techniques  to  the  full  3-D 
compressor.  The  calculations  were  limited  to  extremely  low- 
pressure-ratio  compressors,  due  to  the  perturbation  assumption 
utilized  by  Rae.  Tompkins  et  al . J have  applied  time-dependent 
techniques  to  fully-rotating  compressors. 

More  recently,  Dodge's  method  has  been  extended  to  transonic 
f ully-viscous  solutions  to  stator  systems  having  nearly  irrota- 
tional  vector  diagrams.  It  would  appear  natural  to  extend  this 
method  to  more  general  compressor  cases.  However",  in  the 
presence  of  large  circulation  gradients,  difficulties  confront 
the  direct  application  of  such  an  approach.  The  following  sec- 
tions outline  a numerical  method  for  these  particular  cases. 


Equation  Splitting 


Dodge  closer  ibed  a successful  technique  of  analyzing 
viscous  flow  in  turbomachine  stators.  The  basis  ot  this  tech- 
nique is  substituting  a potential  gradient  plus  a rotational 
component  for  the  velocity  vector. 


W «*  Vij>  + U 


(1) 


substituted  into  the  momentum  Equation 
split  between  potential  and  viscous  terms. 

-VP  + 7 (XV*W)  + V*  In  defW]  l 2) 

e 

where 

P = the  density 
W = velocity  vector 
u * effective  viscosity 
P = static  pressure 

defW  3 V W + (VW)*  where  (VW)*  is  the  transpose  of  VlJ 
\ = second  coefficient  of  viscosity 


When  Equation  (1)  is 
(2)  , the  result  is  a 

DW 

p Pt  = 


Pressure  is  related  to  potential,  as  qiven  by  Equation  (3) 
for  compressible  flow,  leaving  a viscous  Equation  (4). 


where 


V*-V(V$)  = - 

VP 

(3) 

U + pu 

•V(V*) 

- (p*-p)  [V<t»-V  (Vit>)  1 = 

(4) 

( X V • W ) 

+ v. 

Uie  def  W] 

V • V if* 

) 1/y-1 

(5) 

, ( 1_ 

2^I 

Pj  = reference  density 
hj  = reference  stagnation  enthalpy 
V = ratio  of  specific  heats 

Equation  (3)  combined  with  an  equation-of-state  and  the  conti- 
nuity equation  yields  an  equation  for  potential,  which  is 
solved  by  transonic  relaxation.  Dodge  indicates  that  such  a 
splitting,  with  appropriate  boundary  conditions,  results  in  an 
equation  for  which  a potential  vector  satisfying  Equation  (6) 
can  always  be  found. 


‘l  1 


l ' }-h±- 


>-l 


16) 


where 


'1  * l 


>-l 

1 


' h* 

I I 


For  imsos  nivolvuw  littlo  ot  no  spanwise  circulation 
>n'.iiiont  excellent  results  are  obtained  by  such  a method.  For 
example,  Figuie  1 shows  the  suction-surf ace  static-pressure  dis- 
tribution on  a high  Mach  number  stator  recently  tested  in  a 
stage  con?  wurat  ion  by  AFAPL.  Figure  2 shows  the  distribution 
ot  loss  at  the  exit.  Figure  1 shows  calculated  radial-loss  pro- 
files compared  tv'  design  and  two  different  reductions  ot  the 
measured  data.  The  computations  represent  a reasonable  represen- 
t a t i v'n  v't  t he  via  t a , 


However,  when  applying  such  a technivjue  to  c vises  with  largo 
hub-to-shrouvi  circulation  variations,  difficulties  .irise.  For 
example,  consider  axial  flow  in  a cylindric.il  duct  with  a 
constant  axial  velocity  and  no  swirl  in  the  absolute  frame-ot- 
reference  (see  Figure  4).  No  circulation  gradient  exists  in 
the  absolute  reference  frame.  However,  consider  the  rotating 
t rame-ot -ref erence  that  would  occur . The  vorticity  in  this 

is  Jho  whore  w is  the  rotational  frovjuency. 
the  velocity  vector  for  this  case. 


t rame-ot -reference 
Fvpi.it  iv'n  1 7)  viives 


W 


t W 


(7) 


where 


' ■ unit  vector  in  the  v1  direction 
unit,  vector  in  the  a direction 
r * radius 

The  pvifenti.il  component  is  given  by  Kvjuation  (!?),  yielding 
Equations  {'))  . 


„ 4± 

♦ I 34 

3 r 

r §0 

* V 

3r 

+ U 

r 

■>  0. 

1 

— A ♦ U . ■ ..<r 

t vh'  v 

3 


(8) 

(<>a) 
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DESIGN,  CONVENTIONAL  ANALYSIS 

3-D  VISCOUS  CALCULATION 

DERIVED  < ACROSS  THE  BLADE  019  THROTTLE 

DATA  ) THROUGH  THE  BLADE  019  THROTTLE 

USING  STAGE  PEAK  TOTAL  PRESSURE 
BLOCKAGE  DISTRIBUTION  FACTOR  = 1.00 

Figure  3.  Radial  Variation  of  Loss  Coefficient  — 
Comparison  of  Design,  Experiment  and 
3-D  Viscous  Calculation  for  APL 
High-Through-Flow  Stator 
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.'no  approach  is  to  apply  the  inle^  condit  ions  of  Podqe  . I'm 
this.  , the  rotut  tonal  vectoi  II  must  bo  i.lont  ically  aero  t o- 

sul  t inn  in  liquations  (10).  It  will  ho  shown  that  those  equations 
cannot  bo  sa  I i s t t o.l . 


t 1 Ob) 


t l Oo) 


Tito  slat  to  pressure  profile  is  constant,  as  qiven  by  liquation  (III 


l'  I'!.  I 1 - 


7<fr»  V<j>-C2  VY>  '’~1 


W t 1~1 

•’l  (l  - Tfn  ^ 


Hy  intoqratiori  of  liquations  (10b)  and  (10c),  a value  ol  potential 
oan  bo  obtained. 

4*  41  + W 7,  + (u't  ) (rO)  <\  2) 

i)  Z ' 

Piquro  1 shows  a top  view  of  a potential  surface.  Note  that 
the  potential  surface  is  curved  in  the  radial  direction.  Thus, 
the  derivative  of  4 cannot  he  ttoto  in  the  radial  direction. 
Differentiation  of  liquation  (12)  yields  the  same  result. 


-*~t  =*  2utr0 

Thus,  liquation  (10a)  can  not  he  satisfied,  and  the  static 
pressure  is  not  constant  in  either  the  radial  or  0 direction. 
Therefore,  to  obtain  a realistic  solution,  the  rotational 
component  must  he  contained  in  U.  Thus: 

L’()  " uvr  < 


l 

l 5d 


( 1 !>a) 


( l hb) 


( l Sc) 


v. 


equation  (14)  and  equations  (15)  moot  the  requirements, 
but  the  potential  vector  is  severely  misaligned  from  the  velocitv 
vector.  For  a typical  axial  compressor,  the  misalignment  can 
approach  70  degrees.  Such  a consistent  misalignment  will  make* 
the  application  of  normal  transonic  relaxation  methods  impossible 
Thus,  a technique  to  more  closely  align  the  potential  and 
velocity  vectors  must  bo  established.  The  following  is  such  a 
method. 

3-Component  Method 

Replace  the  velocity  vector,  V$,  with  Equation  (16) 

W«^  + Vi(i  + U (16) 


where 


G = initial 

Vd>  = 


guess  at  flow  solution 


A = perturbation  velocity  potential 
u = viscous  perturbation 


A solution  for  potential  is  obtained  by  first  starting  with  the 
momentum  equation  in  a rotating  f rame-of-ref erence  (Equation  17) 


P l (W- V)  W + 2 (to  X W)  - V (C“/2)  J = - VP  + 


(17) 


where 


W 

to 

P 

.2 


= velocity  vector 
= angular  velocity  vector 
= density 


C“  = (tor)2,  where  r2  = x~  + y2 
P = static  pressure 

f.  = shear  and  Reynold's  stress  vector 


Eauation  (17)  is  then  dotted  by  to  yield 
p|w*((W*V)W]  + W»  [2  (to  x W)  ] - 


- W*VP  + wj;  (18) 


The  static  pressure  gradient  may  be  expressed  in  the 
following  form  for  an  ideal  gas  (y  and  R = constant) . 


VP  = a“Vp  + pa 


a<)) 


where 


a = speed  of  sound  (static) 
V = ratio  of  specific  heats 
K = qas  constant 
S = entropy 

Dotting  Equation  (18)  by  W yields 


W*  VP 


(W*Vp)  + pa2  ( Jjj-1  )(W«  VS) 


(20) 


The  continuity  oquation  for  steady,  compressible  flow  is 
given  by’ 


V • ( pW)  = 0 


(21) 


or,  expanding  and  rearranging  Equation  (21)  gives 


W • V p = - p ( V • W ) 


(22) 


Replacing  Equation  (22)  into  Equation  (20)  results  in 


W*VP  = -a2 p ( V • W)  + pa2(^-~)  (W-VS) 


Examining  the  cross-product  term  in  Equation  (18), 

•*  *>  , ► -V 

w x W is  a vector  normal  to  both  W and  m. 

There  fore , 


!3) 


W* (w  x W)  = 0 


(24) 


Substituting  Equations  (23)  and  (24)  into  Equation  (18)  yields 


2 

pjw*  [ (W-V)W  ] - W • V ( ?-)  1 

< r ->  ->  '2t/y-i\- 

= -l-a-p(V*W)  + pa2  ( yr^)(W*  VS) 


(2  5) 


->  -f 

+ vv  • y 


ii 


or  dividing  by  p,  and  rearranging 


W*[(W*V)W]-  W*v(~)  - a~  ( V • W) 

= -a2(^i)  (W-VS)  + i (W 'h 


(26) 


The  velocity  vector,  W,  may  then  be  separated  into  three 
components  as  in  Equation  (16) 


W = G + V<(>  + U 

Substituting  Equation  (27)  into  Equation  (26)  yields 


(27) 


W*[(W*V)(2  + V4>  + U)1  - (2  + V«j>  + U)  *v('2-) 
-a2  (V*  (G  + V4>  + U)  ] 

= -a2^-1)  (W-VS)  + ~ (W*?) 


(28) 


Rearranging  Equation  (28)  , leaving  only  V 4>  terms  on  the  left- 
hand  side,  yields 


W*[(W*V)V4>]  - V<t>*v(C— ) -a2  V • ( V 4> ) 2 

=-W*l(W*V)(6  + U)J  + (6  + u)*v(£-) 


+ a2  (V*  (£  + U)]  -a2 


+ 


i(W-f) 


(29) 


However, 


V*  < V4> ) = V24> 


(30) 
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■and  therefore,  Equation  (29)  finally  becomes 

► 2 
w-  I <W-  V)  J - a2(V2*)  - V*»v(§~) 

2 

* (<5  + i5)  • v(~~)  - w*  ( (w*v)  (d  + tf)  j 

» a2V(2  * U)  - A2(^±)  (W- VS)  * i (W*  t ) (31) 

where  all  terms  on  the  right-hand  side  of  the  equation  are 
presumed  known  from  previous  computations. 

Since  viscous  calculations  are  no^  the  object  of  the  current 
program,  the  appropriate  equation  for  U will  not  be  developed. 
However,  this  development  parallels  that  given  by  Dodge^.  A 
general  outline  of  the  required  numerical  method  is  given  as 
follows : 

o Establish  grid  system 

o Calculate  initial  guess 

o Calculate  relaxation  solution 

These  major  sect  ions  of  the  present  numerical  method  are 
described  below. 

C.rid  System 

The  choice  of  a at  i d system  establishes  many  ot  the  require- 
ments for  the  subsequent  calculations , as  well  as  dictating  the 
information  and  logic  structure  of  the  computer  code.  The 
inlet  circulation  gradient  eonsiderat  ions , discussed  tit  the  pre- 
ceding section,  have  an  impact  on  t It  i s choice.  Consider  the 
grid  system  utilized  by  Dodge1’.  It  is  nonorthogonal  on  a 
succession  of  planes.  Each  plane,  however,  is  orthogonal  to  a 
vector  tield  termed  from  a set  ot  quasi-streamlines.  When 
large  circulation  gradients  are  present,  such  a grid  system  is 
not  possible.  To  demonstrate  this,  consider  the  general  equa- 
tion for  a surface  in  space. 

f (x,y , z)  = C,  ( 32) 


who  re 


is  a constant 


u 


A vector  given  by  Equation  (33) 


dr  = dx  & + dy  9 + dz  z 


lies  on  a surface  when 


3 f . ± 3f  . . . 

v-  dx  + ■;  dy  + dz 

3x  3y  1 3z 


Without  loss  of  qenerality  it  can  be  assumed  that  one  coordinate 
direction  has  a non-zero  derivative  so  that 


- - (§f  ax  * % ay)/  |i 


dr  = dx  x + dy  ? - ( ~ dx  + ~ dy  ) / ~ S 


If  the  surface,  f,  is  to  be  normal  everywhere  to  an  arbitrary 
vector  field,  V,  then 


thus , 


dr • V ' 0 for  all  directions  given  by  dx  and  dy 


ax  vx  + ay  vy  - (Tx 


* H d*)  VM  ‘ 0 


<(vx  H - v2  M ) + dy  ( 


v |1  - v ) = 0. 

y 3z  3y  2 ' 


Since  dx  and  dy  are  arbitrary 


(39a) 
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( 19b) 


" 


an 


V 


3f/3y 

Wit 


* 


it  can  be  shown  that  Equations  (39a)  and  (39b)  hold  only  if 
Equation  (*10)  holds. 

V«Vf  (40) 

Therefore,  there  is  not  a single  surface  everywhere  normal  to  a 
vector  field  unless  that  field  is  irrotational . Likewise,  no 
simple  single  flow  area  exists  describing  any  flow  device,  such 
as  a low-hub- tip-ratio  axial  compressor.  The  view  most  utilized 
by  compressor  designers  is  the  blade- to-blade  section.  Normals 
constructed  across  the  channel  cannot  be  stacked  up  into  anything 
approaching  a surface  orthogonal  to  the  streamwise  direction. 

For  the  modern  low-hub- tip-ratio  comuressor,  this  is  anything 
but  a small  effect.  Large  differences  in  surface  shapes  occur 
between  those  constructed  by  dropping  normals  from  the  mean 
radius  section,  and  those  obtained  by  stacking  blade- to-blade 
normals.  Thus,  one  is  forced  to  consider  a grid  system  and  a 
solution  technique  that  is  completely  nonorthogonal . The  grid 
system  is  selected  by  starting  from  the  mean  section  (hub-to- 
shroud)  on  a user-selected  master  surface,  which  is  also  utilized 
for  radial  equilibrium  calculations.  This  surface,  running 
from  hub  to  shroud,  inlet  to  exit,  but  at  a mean  location 
between  blades,  is  constructed  to  be  two-dimensionally  orthogonal 
(X2-X3  plane,  whore  Xl  is  blade- to-blade , X2  hub- to-shroud , and 
X3  inlet- to-exi t directions).  The  grid  system  is  then  construct- 
ed from  normals  in  tire  blade- to-blade  direction.  This  conforms 
with  the  normal  two-dimensional  practice  in  turbomachinery . 

Blade- to-blade  projections  will  have  the  appearance  (but  not 
the  reality)  of  being  orthogonal.  Hub- to-shroud  views  will 
clearly  be  nonorthogonal  in  cases  where  appreciable  spanwise 
circulation  gradients  exist.  Such  an  example  is  shown  in 
Figure  fa  for  an  axial  inlet  similar  to  that  of  Figure  4. 

Once  the  basic  pattern  of  the  grid  system  has  been  selected, 
the  actual  mechanics  of  the  generation  must  be  developed. 

Dodge  ’ basically  used  a two-step  process. 

(1)  Generate  quasi-streamlines. 

(2)  Carve  the  streamlines  by  projecting  normals 
from  known  points  to  unknown  points. 

Both  processes  consume  considerable  amounts  of  computer  time 
in  the  program  developed  by  Dodge*'.  Some  effort  was  expended 
to  speed  up  this  process,  resulting  in  approximately  a 30  to  1 
improvement.  A brief  discussion  of  this  new  method  follows. 
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A 


MASTER  SURFACE 


A computer  code  can  be  viewed  as  an  information  processor. 
The  code  starts  with  what  is  usually  a very  limited  amount  of 
information  contained  in  the  input  deck.  It  is  then  expanded, 
utilized,  and  ultimately  compacted  back  to  an  amount  small 
enough  for  the  human  mind  to  quickly  assimilate  the  results. 

Input  geometry  typifies  this  process.  The  user  describes  only 
the  boundary  surfaces.  The  computer  fills  out  the  rest  of  the 
grid  system  by  laying  in  a bundle  of  quasi-streamlines. 

However  this  is  accomplished,  the  quasi-streamlines  can  only 
be  a function  of  the  input  data  points.  Equation  (41)  must 
represent  a general  relation  between  a point  on  a quasi-stream- 
line and  the  input  data. 

N 

X.  - £ f..x.  (41) 

1 jti  3 

where 

X^  = coordinate  on  a quasi-streamline 

xj  = input  blade  surface  coordinate 
fji  = influence  coefficient 

The  method  used  by  Dodge*’  to  obtain  quasi-streamlines 
involved  solving  a number  of  2-D  stream-function  problems  to 
first  determine  surfaces  of  revolution  and  then  to  establish 
blade-to-blade  streamlines.  The  current  program  uses  an 
interpolat ive  method.  Surfaces  of  revolution  are  determined  by 
interpolation  between  hub  and  shroud  on  an  equal  area  basis. 
Likewise,  blade-to-blade  surfaces  are  established  by  blade-to- 
blade  interpolation.  It  would  be  expected  that  these  interpola- 
tion factors  could  not  be  independent  for  every  such  surface. 
Thus,  an  efficient  approach  is  to  obtain  a final  set  or  influence 
coefficients  (Equation  41)  by  a convolution  of  hub- to-shroud 
and  blade-to-blade  factors.  Then  the  actual  streamline  points 
can  be  calculated  by  the  lone  vector  multiplication  and  addition 
of  Equation  (41) . Such  an  approach  eliminates  much  of  the 
overhead  that  creeps  into  other  more  fragmented  approaches . 

Once  quasi-streamlines  have  been  obtained,  they  must  be 
cut  into  a grid  system  by  casting  normals  from  selected  stream- 
lines to  neighbors.  The  algorithm  to  accomplish  tills  utilizes 
the  following  steps: 

o Establish  the  order  of  connection  of  quasi- 
streamlines. The  grid  point  on  anv  quasi- 
streamline  is  placed  by  constructing  a normal 
from  a neighboring  quasi-streamline  with  known 
points.  The  order  of  this  connection  in  a 
single  X3-plane  is  shown  in  Figure  7. 
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o Reorder  the  quasi-streamline  definition  to  lino 
up  with  established  grid  connection. 

o Compute  the  fractional  location,  S^,  of  normals 
from  quasi-streamline  points  to  their  connected 
neighbors  and  assure  that  this  array  is  ordered 
(see  Figure  8) . 

o Compute  actual  grid  points  by  establishing  the 

spacing  on  the  master  streamline.  Then  interpo- 
late in  the  table  of  normals  to  find  the  location 
of  the  next  point. 

Such  an  algorithm  minimizes  calculations  by  calculating  normals 
on  grid  lines  that  actually  connect;  calculating  normals  based 
on  the  input  quasi-streamline  points,  which  are  usually  much 
less  in  number  than  grid  points;  and  by  arranging  data  structure 
to  allow  vectorized  loops. 

Initial-Guess  Solution 


The  natural  initial-guess  solution  is  based  on  a 2-D 
solution  on  a hub- to-shroud  surface.  Such  solution  techniques 
were  developed  in  the  early  50' s by  Wu?.  Since  that  time, 
they  have  been  highly  refined  using  both  finite-difference  and 
curvature-of-streamline  techniques.  However,  if  the  mean 
surface  shape  is  specified,  the  governing  stream  function 
undergoes  an  obvious  elliptic-to-hyperbolic  transition  as  the 
Mach  number  passes  through  unity  (in  the  relative  frame) . 
Because  of  this,  solution  techniques  have  concentrated  on 
the  specified  work  formulation,  which  remains  elliptic  until 
the  meridional  component  of  velocity  becomes  supersonic. 
However,  the  work  form  is  inconvenient  when  analyzing  an 
already  existing  geometry,  since  it  is  difficult  to  relate  a 
particular  shape  to  a work  distribution.  For  these  reasons,  a 
transonic  relaxation  solution  to  a Wu-type  stream- function 
equation  was  developed.  An  outline  of  the  technique  follows. 

The  basic  equation  to  develop  a stream- function  solution 
is  the  rotating  version  of  the  Crocco-Vazsonyi  equation  as 
given  by  Owczarek" . 

2(iTxW'  - Wx(VxW)  = - VI  + TVS  + ?/p  (4  2) 


ID 


u)  = rotational  vector 
I = rothalpy 
T = static  temperature 
S = entropy 

? = viscous  stress  vector 
p = density 
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Wu  concentrated  on  cylindrical  coordinates.  To  be  useful 
for  all  configurations,  an  equivalent  form  in  arbitrary  curvi- 
linear coordinates  is  required.  Only  one  component  of  Equation 
(43)  is  actually  utilized;  in  this  case,  the  component  in  the 
X2  direction.  It  is  assumed  that  flow  only  exists  parallel  to 
the  master  surface,  yieldinq  W j 0. 


W2  8h2l 

h2)  V 3 ;ix3 


T 3S 1 31  . ■'  2 

il,  h.,  p 


(4  ]) 


A stream  function  can  be  defined  that  automatically  solves 
continui ty . 


1 

h.„ 


t,DW3 


(44) 


where  n is  a factor  to  take  into  account  blade-tc-bl ade  stream- 
sheet  convergence. 


Prom  these  equations,  one  can  follow  Wu's  development  directly, 
yielding  a stream-function  equation  qiven  by  Equation  (46) 


2 ■?  w w ')  o 

2 w2  1 n1  -W2Wj  y-y  2 2 l 32y 

(a  ' 2)  h|  ax2  " l‘2TS  h2  ax2 


+ + M.  ?y.  = o 

h,  DX2  h3  DX3 


where 


2 2 
N = (a  -W  ) 


+ 2W,u,  - 


M = (a2-W2) 


(l 

T2 

3S  1 
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staqnation  enthalpy  at  the  streamline  inlet 


Oo 


(46) 


absolute  preswirl 


. H7 


The  solution  procedure  ls  essential ly  similar  to  that  for 
potential  equations.  The  coefficients  are  computed  from 
Equation  (46) . The  result  is  then  a coefficient  matrix  whose 
region  of  influence  is  controlled  by  the  local  Mach  number  with 
boundary  conditions  which  specify  4'  on  all  boundaries.  The 
stream  function  is  then  differentiated  to  form  a pW  product 
which  can  be  two-valued  in  W.  Both  supersonic  and  subsonic 
solutions  are  obtained  by  nondimensionalizing  the  pW  curves  to 
the  form  given  by  Equation  (47). 


G 


,)W 

T’TT 


TT  - 


-21-  (i 
/r=l  V 


vL/y-l 


(47) 


where 


W*  = U/>T?r 


p"  = stagnation  density  in  the  relative  frame 
h"  = stagnation  enthaloy  in  the  relative  frame 
a"  = /(V-l)hH 


Both  the  subsonic  and  supersonic  velocities  are  obtained 
from  approximation  to  an  inverse  of  Equation  (47),  as  given  by 
Equation  (48)  for  subsonic  flow  and  (49)  for  supersonic  flow. 


W*  = W* 
1 m 


1 - (1-G/G 


m 


W*  = /2- 


W* 

m 


1 — (,/G  t W* 


m 


m 


(48) 

(49) 


where 


_ Y + l 

, ,2rFTT 


( I — 

m 
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W = 
m 


2 (Y-l 


V Ttr 


G and  W represent  the  maximum  value  of  flow  and  the 
m m 

velocity  at  this  maximum,  respectively. 
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One  must  still  decide  which  solution  to  select.  This  is 
accomplished  based  on  a direct  analoqy  to  potential  solution. 
A potential-like  parameter  is  computed  along  streamlines  from 
inlet  to  exit  by  computinq  the  following: 

Exit 

r { T ) = / W h3dx3  ( 

Inlet 


The  following  conditions  can  be  applied  to  W,  depending  on 
the  case. 


Case 

Subsonic 

Throughout 

Supersonic  In 
Subsonic  out 


Supersonic  Tip 
Subsonic  Hub 


Subsonic  In 
Supersonic  Out 

Supersonic  In 
Supersonic  Out 


1'  (40 


Irrelevant — always  take  W^* 


React  a desired  value  of 
r(Y),  make  a transition 
from  W2*  to  Wj*  at  the 
point  yielding  an  exit 
integral  of  I’ (S')  equal  to 
that  desired. 

On  hub  use  subsonic  value 
throughout.  An  exit- 
circulation  spanwise  I’ 
variation  is  then  specified 

Make  transition  at  minimum 
area  location 

Use  only  supersonic  values 


With  the  specification  of  these  conditions,  a solution  to  the 
flow  on  the  mean  stream  surface  is  provided.  For  cases  with 
supersonic  in  and  subsonic  out,  the  actual  shock  location  is 
arbitrarily  dependent  on  the  choice  of  the  circulation  parameter. 
In  these  cases,  the  initial  guess  may  be  too  realistic.  To 
illustrate,  consider  the  equivalent  of  Equation  (31)  for  flow 
in  a 1-D  duct  with  area  variations. 
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(w2-a2)  ^4  ~ 

dx~ 


2 d$  1 dA ' 
1 ilx  /Tr  dx 


- d (I 


2 1 dA' 

° A1'*  3T" 


- (W-a~) 


where 

l dA1  = 1 dA  _ / >-l \ d£ 
A7"  dlT-  A dx  V > R / dx 


A is  the  cross-sectional  area 

For  this  case,  entropy  changes  are  indistinguishable  from 
area  variations.  Considering  the  diverging  portion  of  a 
converging-diverging  nozzle  at  a high  back  pressure,  a shock 
stands  in  the  nozzle.  Flow  upstream  and  downstream  of  this 
shock  approximately  obeys  the  1-D  equations.  Across  the  shock 
the  Rankine-Hugoniot  relations  apply.  The  actual  location  of 
the  shock  is  fixed  by  establishing  the  precise  points  where  the 
Rank ine-Hugoniot  relation  is  in  balance  with  upstream  and 
downstream  solutions.  If  the  guess  is  very  good,  such  as 
solving  isentropic  relations  both  ahead  of  and  behind  a discon- 
tinuity arbitrarily  located,  then  any  adjustment  in  shock 
position  must  be  accomplished  by  recognizing  the  error  at  the 
shock . This  error  must  then  propagate  to  the  whole  flow  field, 
as  indicated  schematically  in  Figure  9.  There  are  several 
problems  in  obtaining  such  a solution.  The  first  of  these  is 
requiring  an  isentropic  rise,  .is  well  as  tall,  in  the  potential 
solution.  The  star-switching  methods  in  use  discriminate 
against  such  potential  rises.  Additionally,  the  error  is  only 
at  the  guessed  discontinuity  and  it  must  propagate  out  to  the 
test  of  the  field.  Figure  10  shows  the  results  that  actually 
occur  for  such  a guess  (1st  calculation).  Hxcept  in  the  region 
ot  the  discontinuity,  little  potential  correction  occurs.  A 
much  better  choice  is  to  connect  input  to  output  by  a straight 
line  yielding  continuous  flows  and  relatively  large  differences 
between  guessed  and  final  solutions  everywhere.  This  approach 
is  used  for  the  full  three-dimensional  solutions  described  in 
Section  II L. 


The  mean  surface  solution  is  projected  to  the  surfaces  by 
utilizing  Equation  (52). 


% 


u 


Figure  9.  The  Effect  of 


a Good  Initial  Guess 


1 component  ot  velocity 


grid  spacing  in  the  1 direction  on  the 
master  surface 

h.  grid  spacing  in  the  l direction 
C ^ l component  of  guessed  velocity 


Mote  that  such  a projection  assures  that  the  solution  preserves 
blado-to-blado  circulation.  Cl  > is  simply  made  equal  to  the 
mean  surface  solution  value,  and  C.,  is  assumed  aero.  Ill 
losses  are  incorporated  in  the  guessed  solutions  through  an 
input  loss  and  blade  boundary-layer  blockage  distribution. 


To  this  component  one  must  add  a rotat ional  component  asso- 
ciated with  radius  change.  This  can  be  derived  from  the  rotat mg 
reference  frame  equivalent  of  the  Crocco-Vaz  sonv  i , equation  (4.’). 
The  result  for  the  t'  component  of  equation  (-Id)  is  given  below. 


31  v 3W  0 

JT  - (W0+  uir)  sina  + r ~- 


where 

W is  the  meridonal  velocity 
ni 

W ^ is  the  tangential  velocity 

ui  is  the  rotational  frequency 

a is  the  streamline  slope 


When  o is  zero  the  above  equation  is  irrotational . However,  when 
o is  non-zero  a rotational  component  ex i s£s  that  is  essentially 
geometric  and  must  be  included  into  (1  or  U since  cannot  exhi- 
bit rotational  behavior.  This  component  is  approximated  by  a U 
component  in  the  streamwise  direction  as  given  below. 


where 


U 


3 


or 

cos  P 


sina 


(t» 


CL 


) 


t>  is  the  grid  system  angle  in  the  streamwise  direction 

i a >.  - 1 , rd  t) . 

given  by  = tan  ) 

0 is  the  value  of  0 on  the  master  surface  across 
the  blade  passage 


This  component  is  wholly  geometric,  zero  on  the  master  surface, 
and  is  simply  added  to  the  results  ot  liquation  (T>J). 


Potential  Solution 


The  potential  solution  proceeds  in  a similar  fashion  to 
that  described  by  Dodged.  However,  the  use  of  a completely  non- 
orthogonal  grid  seriously  complicates  the  geometry  considera- 
tions . 


Equation  (31)  represents  a general  vector  form  of  the 
partial  differential  equation  for  potential,  which  may  now  be 
expressed  in  Cartesian  coordinates  as  follows. 
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Equation  (53)  represents  the  partial  differential 
equation  for  perturbation  potential  in  Cartesian  coordinates. 

The  equation  is  elliptic  for  subsonic  conditions  and  hyperbolic 
for  supersonic  conditions.  Thus,  it  is  adaptable  to  solution 
by  relaxation  techniques. 

Boundary  Conditions  on  Perturbation  Potential  - At  the  inlet 
and  outlet  of  the  blade  passage,  the  perturbation  potential,  <p,  is 
assumed  to  be  equal  to  zero.  This  is  a valid  assumption  when 
the  inlet  and  outlet  are  chosen  far  enough  from  the  blades  that 
the  radial  equilibrium  solution  provides  a satisfactory  repre- 
sentation of  the  flow.  On  solid  surfaces,  and  bounding  stream- 
lines, the  boundary  condition  on  $ is 

V <+>  • n = 0 ( 54 ) 

or,  expanding  in  Cartesian  coordinates 

3 <t>  , 9 <J>  . _ 3 4>  _ 

n 7^-  + n -r-A  + n -r-=r)  (55) 

x 3x  y 3y  z 3z 

Car tesian-to-Nonorthogonal  Coordinate  Transformation  - 
Equation  ( 5 3)  Is  expressed  in  terms  of  Cartesian  coordinates ; 
however,  the  grid  system  generated  within  the  program  is  non- 
orthogonal.  Thus,  in  order  to  formulate  difference  schemes  in 
the  nonorthogonal  grid,  the  Cartesian  derivatives  in  Equation 
(53)  must  be  recast  in  terms  of  derivatives  with  respect  to  non- 
orthogonal coordinates.  A transformation  matrix  to  accomplish 
this  was  developed  as  follows. 

Consider  a scalar  property  A.  The  differential  of  A with 
respect  to  Cartesian  coordinates  is: 


..  3 A . 3A  . 

dA  = -s—  dx  + -r—  dy  + 

3x  9y 


, 3A  , 

dz  = 2^  w aX^ 


where 


Equation  (56)  provides  the  basis  for  expressing  derivatives 
in  the  nonorthogonal  coordinate  directions  (xj_,  X2,  X3)  . The 
partial  derivative  of  A,  with  respect  to  the  nonorthogonal 
direction  xi,  may  then  be  formed  from  Equation  (56) : 


3 ax. 

3 A 3 A j _ . t . 

dxi  3X7  3x7  for  1 “ lf  2 9 3 
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Second  derivatives  may  be  treated  similarly,  with  the 
differential  written  in  Cartesian  coordinates  as: 


0 I#-)  ■ £ hr  (Hr)  d\  1 * l>  2' 

1 k=l  k l 


Then,  the  partial  derivative  or  A—  , with  respect  to  the  non- 
orthogonal  direction  x,,  becomes  ' 1 


(3A 
lx,  Ux 


V-  3 / 3 A \ k ^ 1 o l 

£ Tv-  a ft'  777"  for  1 = 3'  3'  3 


, \ lx  . ' ' IX.  \ lx.  / lx  0 

t i k=l  kit 


1*1/  2,  3 


Substituting  Equation  (ri7)  into  Equation  (59)  yields  the 
following : 


. 3 , r 3 , IX.  a IX, 

*L.)  = v 3 y 3 i 

x.  1 IX,  A.  IX.  lx.  lx 

i k=l  k L ] = 1 3 l J 


lx  ' lx 
x» 


or,  expanding  Equation  (60) 


ix  , \ ix. > .a,  \ ix.  ix.  ' 7x7 

t i k=l  j=l  k 3 it 


3 3 v , IX.  v IX. 

+ y y / 3 A_  3 ( i k 

k=l  3=1  3 k i t 


for  i = 1,  2,  3 
1 “ 1/  2,  3 


In  Equation  (61),  the  geometry  derivative  -y^—  ’HT  must  be 

‘ k i 

recast  in  terms  of  nonorthogonal  derivatives.  This  conversion 
is  performed  by  rewriting  the  second  term  of  Equation  (61)  as: 


and  then  employing  Equation  (57)  to  replace  the  term  in  brackets: 


Then,  the  final  form  of  Equation  (61)  is 
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A 12  x 12  transformation  matrix  may  be  built  to  represent 
Equation  ( 62) : 


3X 


(9  x 9) 


— 0 — 


(3  x 9) 


i_  (-11) 

'X  £ ^x“ 


(9  x 3) 


(3x3) 


Defining  the  12-element  derivative  vectors  as 


Equation  (62)  is  represented  as 
= TD 


I 

Then,  by  taking  the  inverse  of  T,  the  derivatives  with  respect  to 
Cartesian  coordinates  may  be  expressed  in  terms  of  derivatives 
with  respect  to  nonor thogonal  coordinates: 


Equation  (61)  thus  provides  the  necessary  transformation 
I expression  for  the  conversion  of  derivatives  in  Equations  (03) 

and  (55)  to  nor.-' r nosonal  coordinates. 

I 

I Numerical  Motive* ■-  for  Solution  of  Potenti.il  Equation  - The 

partial  dif lerential  equation  governing  perturbation  potential 
[Equation  (53)]  is  recast  as  a linear  difference  equation  in  non- 
orthogonal  coordinates.  The  difference  technique  is  described 
in  Appendix  A.  The  resulting  linear  system  may  then  be  solved 
by  relaxation  techniques. 

The  difference  star  used  employs  40  nodes,  and  is  capable 
of  treating  derivatives  at  all  interior  ard  boundary  nodes.  It 
produces  a centered  difference  scheme  in  the  streamwise  direction 
for  subsonic  flow,  and  a mixed-hyperbolic  and  implicit  scheme 
in  the  presence  of  supersonic  flow.  The  degree  of  mixing  is 
determined  as  a function  of  local  Mach  number,  based  upon  a 
criterion  developed  from  a stability  analysis  of  the  basic 
relaxation  difference  equation.  As  Mich  number  increases,  the 
variable  supersonic  difference  star  becomes  weighted  to  more 
closely  reflect  the  region  of  influence  of  the  characteristic 
lines,  eventually  becoming  totally  hyperbolic. 

Program  Architecture  - In  designing  this  program,  a primary 
goal  was  to  maintain  the  capability  of  handling  problems  to  a 
maximum  size  of  10,000  nodes.  In  order  to  accomplish  this  on 
the  CY174,  it  became  necessary  to  introduce  sequential  files 
within  the  relaxation  module  to  handle  large  data  blocks,  such 
as  grid- transformation  arrays  and  coefficient  arrays.  In 
addition,  smaller  data  blocks  were  stored  on  random  files 
managed  by  the  scratch  information  processor.  A general  logic 
diagram  of  the  potential  equation  solution  module  is  given  in 
Figure  11. 

The  relaxation  iteration  package  utilized  in  this  program 
consists  of  an  optimized  assembly  language  routine,  which  per- 
forms ordinary  over-relaxation  and  permits  the  entire  problem 
to  be  held  within  core  memory  during  iteration.  The  assembly 
language  routine  is  more  than  three  times  faster  than  the 
Fortran  routine  that  it  replaces. 
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SECTION  III 
RESULTS 

The  USAF  has  had  a program  to  develop  high-pressure-ratio, 
high-diffusion  axial  compressors  for  several  years.  In 
particular,  the  program  of  interest  involved  a rotor  designed 
for  a 3.0  to  1.0  pressure  ratio  in  a 1600  ft/sec  axial  stage. 

Initial  tests  of  this  device  indicated  that  deviation  angles 
were  excessive  without  a second  partial-flow  splitter  near  the 
exit.  A center-section  from  this  rotor  was  tested  with  a 
splitter  and  the  results  were  reported  by  Holtman  et  al-  . 

Dodge1^  later  studied  this  cascade  utilizing  two-dimensional 
transonic  relaxation  methods.  The  rotor  was  designed  by 
Wennerstrom  and  Frost11,  and  the  test  results  were  reported  by 
Wennerstrom  et  al-*-2.  Although  the  computer  program  developed 
herein  is  expected  to  have  much  wider  application  than  these 
devices,  the  motivation  for  its  development  rests  in  certain 
inabilities  to  design  such  high  work  stages. 

Although  the  splitter  had  substantial  benefits,  its  cascade 
test  indicated  that  rather  severe  diffusions  occur.  However, 
questions  existed  as  to  the  applicability  of  the  cascade  test. 

The  extreme  convergence  of  the  stage  required  contracting 
relatively  low-aspect-ratio  end  walls  applied  to  a supersonic 
cascade.  This  conf iguration  introduces  a complex  interaction 
between  curvatures  simultaneously  occurring  in  two  directions. 

The  intent  of  this  program  was  to  explore  these  questions  by 
analysis.  The  results  of  this  study  are  described  below, 
beginning  with  a conventional  transonic  cascade,  then  high- 
turning  cascades,  and  finally  a rotor. 

VKI  Low-Camber  Cascade 

In  1973,  as  part  of  a von  Karman  Institute  short  course, 

Breugelmans  and  Starken12  reported  data  on  a low-camber  DC  A 
airfoil.  As  a test  case,  a two-bladed  cascade  configuration 
was  calculated  at  an  inlet  Mach  number  of  1.409.  The  grid 
system  along  the  centerline  is  shown  in  Figure  12.  Note  that 
the  two  side  walls  are  formed  by  the  suction  and  pressure  sur- 
face of  the  cascade.  The  resultant  solution  is  shown  in 
Figure  13.  The  center  blade  and  outside  wall  show  very  similar 
pressure  distributions.  However,  the  pressure  surface  of  the 
center  blade  indicates  some  of  the  previously  unpredictable 
acceleration  near  the  leading  edge.  Two  data  cases  bracketing 
the  calculated  static  pressure  ratio  are  compared  to  the  center 
blade  loading  in  Figure  14. 

The  one  very  low  static  pressure  point  on  the  suction 
surface  just  ahead  of  the  trailing  edge  is  considered  spurious, 
since  it  does  not  appear  on  any  other  sections.  Agreement  is 
as  good  as  the  comparisons  between  time-dependent  and  non- 
orthogonal  relaxation  methods, presented  by  Dodge1  previously, 
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Figure  12.  Centerline  Grid  System  VKI  Low- 
Camber  DC A Cascade,  Mi  = 1.409, 
32  = 47.0  Degrees 
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Figure  13.  Static  Pressure  Ratios  Calculated 
on  Two-Blade  VKI  9.5-Degree 
Camber  Cascade,  = 49.0  Degrees 


r 


with  the  leading  edqe  regions  reproduced  better.  Principal 
regions  of  disagreement  are  in  the  exit  region — a region 
dominated  by  viscous  effects. 

AR1.  Cascade 

The  ARL  cascade  has  been  tested  in  two  forms.  The  first 
was  without  a splitter  tested  by  Holtman  et  al*4,  and  the  second 
was  with  splitters  by  Holtman  et  al0.  A detailed  study  utilizing 
a two-dimensional  transonic  relaxation  method,  described  by 
Dodge 1 , has  already  been  reported.  Results  agreed  well  with  the 
data.  However,  viscous  effects  had  to  be  modeled  by  reducing 
the  normal  stream  tube  convergence  (b-width) . This  is  required 
since  the  calculation  technique  is  inviscid.  The  result  was  that 
the  real  performance  map  of  the  cascade  could  not  be  matched 
because  it  was  throttled.  The  effect  of  this  deficiency  is 
shown  in  Figure  15,  plotting  axial-velocity/density  ratio  versus 
exit  Mach  number.  To  close  exit  loading  properly,  the  correct 
axial-velocitv/density  ratio  must  be  selected.  However*,  to 
reasonably  match  surface  loadings,  the  exit  Mach  number  must  be 
duplicated.  These  two  quantities  are  far  apart  for  this  cascade 
when  an  inviscid  technique  is  utilized.  The  application  of  the 
current  method  alleviates  this  problem.  Results  are  shown  in 
Figure  16  for  surface  static-pressure  distribution. 


Two  significant  differences  between  experiment  and  calcula- 
tion appear.  A compression  on  the  suction  surface  followed  by 
an  acceleration  with  a near- trail ing-edge  compression  is  pre- 
dicted but  not  evidenced  in  the  data,  or  in  any  of  the  calcula- 
tions of  Dodge*.  The  origin  of  this  phenomenon  is  not  yet  clear. 
However,  it  is  associated  with  the  high  loading  of  the  split- 
tered  cascade,  since  it  disappears  with  the  lower-loaded  split- 
tered  cascade.  A Prandtl-Meyer  expansion  with  20  degrees  of 
turning  would  produce  such  an  acceleration,  while  the  turning  on 
the  suction  surface  over  the  same  range  is  35  degrees.  However, 
it  would  be  unlikely  to  occur  in  a test  because  of  the  separa- 
tion after  the  first  compression,  and  it  may  not  have  appeared 
in  earlier  calculations  because  of  the  lack  of  three- 
dimensional  i ty  of  previous  tests.  A complete  contour  map  of  the 
suction  surface  is  shown  in  Figure  17. 


The  other  discrepancy  is  well  documented.  A similar  but 
smaller  accelerative  peak  occurs  on  the  pressure  surface  near 
the  leading  edge.  Dodge*®  extensively  studied  this  phenomenon , 
and  came  to  the  conclusion  that  it  was  due  to  pressure  surface 
curvature,  alternating  between  expansive  at  the  very  leading 
edge  and  compressive  further  down.  Figure  18,  reproduced  from 
Dodge*®,  shows  that  the  pressure  surface  crosses  over  in  the 
vicinity  of  20  to  30  percent  of  chord.  This  effect  persists 
through  all  calculations  by  all  techniques,  and  thus  must  be 
considered  to  be  real.  However,  it  does  not  appear  in  data  due 
to  rapid  boundary-layer  thickening. 
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Figure  17.  Suction  Surface  Contours  of  Constant  Static 
Pressure  Ratio  for  ARL  Cascade  Without 
Splitter,  P2/P1  =1.88 
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Figure  19  shows  the  blade  shape  of  the  splittered  cascade. 
Like  the  unsplittered  cascade,  it  was  tested  as  a two-dimensional 
blade  with  a sharply  converging  end  wall.  Figure  20  shows  cal- 
culated surface  static-pressure  distributions.  Figure  21  shows 
the  data  for  the  equivalent  case.  Comparison  between  the  data 
and  calculations  are  shown  in  Figures  22  through  24.  Figures  25 
through  28  show  the  three-dimensional  calculated  surface  static 
pressures.  The  splittered  cascade  is  more  three-dimensional 
than  the  unsplittered  cascade,  as  shown  by  comparing  Figure  25 
to  Figure  17.  The  assumed  centerline  loss  distribution  is 
shown  m Figure  29.  This  is  a highly  three-dimensional  case 
without  the  traditional  cascade  practice  of  a large  uniform 
center  span.  It  would  be  nearly  impossible  to  assess  what  the 
effective  two-dimensional  b-width  variation,  so  familiar  to 
conventional  compressor  design,  would  be  for  such  a case. 

However,  by  comparison  to  the  equivalent-  rotor  section,  the 
validity  of  the  cascade  approach  would  be  verified.  Agreement 
with  the  data  is  excellent.  The  pressure  surface  also  shows 
the  acceleration  in  the  vicinity  of  the  splitter  leading  edge. 

The  previous  study  by  Dodge'*' ^ was  only  successful  for  b- 
width  distributions  near  the  metal  (those  without  end-wall  bound- 
ary- layer  corrections)  rate  of  convergence.  On  a one-dimensional 
area  basis,  it  could  be  demonstrated  that  the  lower  passage  was 
choked.  When  the  calculation  is  performed  on  a fully  three- 
dimensional  basis,  as  demonstrated  herein,  solutions  are  possi- 
ble. This  is  apparently  due  to  the  combination  of  using  two 
three-dimensional  areas,  flow  readjustment  from  lower  to  upper 
spJ ittered  passage,  and  the  specified  loss  instead  of  b-width 
reduction . 


ARL  Rotor 


Calculations  were  made  on  the  ARL  rotor  described  by 
Wennerstrom  and  Frost^.  Tests  on  the  rotor  were  inhibited  by  a 
stator  choke,  which  caused  the  rotor  to  run  at  very  low  flow. 

The  effect  was  most  pronounced  at  the  tip,  where  deviation 
angles  exceeded  30  degrees.  Total  pressure  loss  also  approached 
50  percent  at  the  tip.  Boundary  conditions  for  the  calculation 
were  taken  from  Wennerstrom  et  al^-,  for  the  through- the-blade 
match  to  the  100-percent  speed  design  point.  Projections  of  the 
grid  systems  at  the  hub,  mean,  and  tip  are  shown  in  Figures  30 
through  32.  Calculated  surface  static  pressures  are  shown  in 
Figures  33  through  36.  At  high  flow  angles,  flow  conditions  are 
extremely  sensitive  to  inlet  and  exit  stream-function  boundary 
conditions.  Thus,  the  exit  conditions  were  not  precisely 
matched.  Inlet  Mach  numbers  are  generally  lower  and  flow  angles 
higher  than  design  because  of  the  very  low  flow  condition. 
Conditions  are  not  as  extreme  as  might  be  expected,  since  the 
tip  passes  very  little  flow.  Loss  distribution  through  the 
blade  was  set  to  linear  for  lack  of  any  better  alternative. 

This  leads  to  relatively  high  pressures  at  the  tip.  These  are 
probably  not  realized  in  the  experiment  due  to  very  large 
near- leading-edge  total-pressure  losses.  Blade  sections 
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Figure  24.  Splitter  Blade  Static  Pressure 
Distribution,  P0/P^  = 1.88 
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generally  possess  an  expansion  near  the  leading  edge,  raising 
suction  surface  Mach  numbers  to  design  or  above.  Direct  compari 
son  was  made  to  the  mid-line  section,  since  it  was  the  one 
modeled  in  cascade  testing.  These  surface  loadings  arc  shown 
in  Figure  37.  Note  that  suction-surf ace  static  pressures  are  as 
low  as  the  cascade,  despite  the  lower  Mach  number  upstream  (1.34 
rotor  to  1.46  cascade). 

The  exit  static  pressure  is  slightly  higher  than  reported 
for  the  experiment.  Rather  than  rerun  the  calculation  of  the 
rotor,  a cascade  calculation  already  available  at  the  higher 
static  pressure  was  used.  The  results  are  shown  in  Figure  38. 
The  comparison  between  the  rotor  and  experiment  is  shown  for  the 
main  blades  in  Figure  39,  and  for  the  splitter  in  Figure  40. 

The  agreement  is  excellent,  providing  strong  evidence  that  the 
cascade  was  a good  simulation  to  the  center  section  of  the  rotor 
This  fact  was  already  somewhat  experimentally  confirmed  by 
comparing  rotor  cascade  loss  coefficients. 
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Figure  40  Comparison  of  Splitter  Blade  Static  Pressure 
Between  Cascade  and  Rotor  Mean  Section 
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SECTION  IV 

CONCLUSIONS  AND  RECOMMENDATIONS 

The  most  significant  conclusion  from  this  study  is  that  a 
technique  has  been  developed,  coded,  and  to  limited  degree 
verified  offering  substantial  insight  into  the  three-dimensional 
flow  fields  internal  to  compressors.  Moreover,  this  computer 
code  is  a practical  design  tool.  Run  times  range  from  10  to  12 
minutes  on  the  CDC  174  (1.0  to  1.2  minutes  on  the  CDC  7600). 

Input  is  arranged  so  that  almost  any  conceivable  compressor 
configuration  can  be  handled,  including  multiple  splitters  of 
varying  shape  and  extent. 

Comparison  between  the  calculation  and  the  splittered 
cascade  tested  by  ARL  is  excellent.  When  the  calculation  is 
applied  to  the  rotor,  the  section  that  corresponds  to  the 
cascade  in  shape  and  overall  construction  has  almost  precisely 
the  same  predicted  static  pressure  distribution  as  the  cascade 
prediction.  This  suggests  that  the  cascade  is  reasonably 
representative  of  the  rotor  section.  A word  of  caution.  This 
conclusion  only  applies  to  this  particular  section  through  this 
particular  rotor.  This  section  was  carefully  selected  by  the 
designers  to  have  minimum  three-dimensional  effects.  Other 
sections  and  other  rotors  need  to  be  examined  on  a case  by  case 
basis,  which  is  now  possible  with  this  new  numerical  method. 

The  close  agreement  between  computation  and  experiment 
tends  to  support  the  concept  of  a supersonic  inlet  cascade  with 
sharply  coverging  end  walls.  The  end  wall  has  a powerful  effect 
on  the  flow  distribution;  however,  the  centerline  results  indi- 
cate reasonable  simulation  of  stream-tube  convergence  is  present 
in  the  rotor.  Again,  care  must  be  exercised  against  over 
generalizing  this  conclusion.  Each  new  case  needs  to  be  exami- 
ned on  its  own  merits. 

Several  years  of  intensive  development  of  massive  numerical 
computer  simulations  to  internal  flow  devices  has  shown  that 
performance  may  be  strongly  affected  by  choices  of  the  tradi- 
tional overall  parameters;  such  as  solidity,  aspect  ratio, 
diffusion  ratio,  etc.  However,  after  these  parameters  are 
optimized  to  a particular  application,  the  designer  still  faces 
many  decisions  on  the  final  blade  shape.  These  choices  are  not 
trivial.  They  can  and  do  affect  performance  substantially;  in 
some  cases,  a great  deal  more  than  do  the  traditional  parameters. 

Without  internal-flow  visibility,  total  optimization  is  not 
possible.  In  the  case  of  turbine  stators,  where  significant 
differences  in  surface  static  pressure  are  not  as  likely  to 
occur,  internal-flow  visibility  must  take  the  form  of  completely 
viscous  solutions.  However,  a brief  perusal  of  any  high  Mach 
number  compressor,  such  as  the  one  shown  herein,  indicates 
significant  regions  of  unwanted  diffusion.  Much  better  per- 
formance could  be  obtained  by  eliminating  such  diffusion.  This 
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study  presents  the  method  to  predict  and  thereby  eliminate  un- 
necessary diffusion.  The  ultimate  test  should  now  be  experi- 
mental. An  extensive  optimization  of  a stage  with  already  good 
performance  will  indicate  what  the  final  payoff  in  design  of 
this  numerical  method  will  be. 
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APPENDIX  A 


RELAXATION  METHOD 


The  major  extension  of  the  current  program  over  that 
reported  earlier  by  Dodge6  is  in  the  relaxation  area.  Since  the 
discovery  of  transonic  relaxation  techniques  by  Murman  and 
Cole15,  the  technique  for  solving  Equation  (31)  has  been  avail- 
able. Since  that  time,  much  has  been  learned  about  this  tech- 
nique. No  attempt  to  discuss  relaxation  processes  in  detail  is 
contained  herein.  Rather,  an  explanation  of  what  is  in  the 
current  program  and  a brief  discussion  of  why,  follows. 


In  its  simplest  form,  Murman ' s method  provides 
approach  for  solving  Equation  (A-l). 


u-M2)  . r±  , a i!*  + . o. 

3s  3n^  3S  3n 


a consistent 


(A-l) 


2 

When  M is  less  than  unity,  Equation  (A-l)  is  elliptic.  Ordinary 
centered  differences  result  in  a relaxation  process  essentially 
equivalent  to  that  used  in  the  incompressible  program.  When 
M is  greater  than  one,  the  equation  becomes  hyperbolic,  and  the 
region  of  influence  changes  to  that  contained  by  a set  of  char- 
acteristics. To  model  this,  differences  must  be  switched  to 
backward  at  the  sonic  point.  For  supersonic  flow,  the  difference 
Equation  (A-2)  was  given  by  Murman  and  Cole15. 
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where  s = s + iAs 

o 

n = n + jAn 
o J 

If  the  difference  solution  for  potential  is  assumed  to  be  a 
continuous  function,  it  can  be  expanded  by  a Taylor's  series 
about  >.*’1,  j.  When  this  is  done  and  substituted  into  Equation 
(A-2),  Equation  (A-3)  results. 
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Note  that  the  equation  is  consistent.  The  lowest  order  trunca- 
tion is  then  given  by  Equation  (A-4) . 


T = - (M  — 1 ) 


ill 

3s3 
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(A-4) 


In  the  conventional  sense,  this  difference  is  first  order  in  s. 
Note,  however,  that  the  magnitude  of  truncation  error  increases 
with  Mach  number.  At  unity,  when  the  difference  equation  and 
the  differential  equation  match  regions  of  influence  exactly, 
the  system  is  second  order. 

Dodge3-  suggested  a nonorthogonal  method  that  was  second 
order,  and  matched  precisely  regions  of  influence.  In  a fixed 
grid  this  method  can  be  applied  in  an  approximate  fashion. 
Consider  the  difference  Equation  (A-5) . 
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The  equivalent  partial  differential  equation  to  Equation  (A—  5 ) 
is  given  by  Equation  (A-6) . 
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where  m + r = k 

The  first  few  terms  of  this  series  are  given  by  Equation  (A-7) 
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If  Equation  (A-8)  holds,  then  Equation  (A-7),  and  thereby 
Equation  (A-5) , are  second-order  consistent  approximations  to 
Equation  (A-l) . 


(A-8) 


At  this  point,  the  characteristics  of  the  differential  equation 
and  difference  equation  match  precisely,  as  indicated  below. 


However,  this  is  only  a consistent  approximation  when  Equation 
(A-8)  applies.  An  approximation  that  is  second-order  when  M is 
unity  and  when  Equation  (A-8)  holds,  and  first-order  in  between, 
is  given  by  Equation  (A-9) . 


*4 


i-1,  j+1 


+ *4 


i-2,  j 


+ 


♦i,j+i  ' 2h,j  + 


(A-9) 


Substitution  of  Equation  (A-9)  into  Equation  (A-l)  results  in  a 
consistent  approximation. 


In  the  full  3-D  case,  first  derivatives  and  mixed  second 
partials  must  also  be  considered.  Instead  of  using  centered 
differences,  the  following  are  used: 
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where  A 1 or  0 

A is  unity  when  the  eoetticiont  nl  (ho  (liaijon.il  iornis 
divided  by  the  eoetticiont  of  the  mixed  partial  makes  an 
opposite  contribution  to  the  diagonal.  Thus,  the  ditti'renee 
contributes  to  the  diagonal  dominance  ot  the  result inq  matrix. 
A similar  lorm  is  used  for  first  derivatives. 


Alter  trying  a number  of  different  block  relaxation  and 
alternatin'!  direction  methods,  ordinary  over-relaxation  was 
t inally  selected  for  the  actual  relaxation  process.  This  is 
implemented  in  a special  assembly  language  coded  routine  that  i: 
optimized  for  storage  and  speed.  The  routine  is  more  than  throe 
times  faster  than  the  Fortran  routine  it  replaces. 


GLOSSARY  OF  TLRMS 


A Area 

a Speed  of  sound 

C Local  wheel  speed 

G Guessed  velocity 

h Lnthalpy 

h'  Total  enthalpy 

I Rothalpy 

P Static  pressure 

P'  Total  pressure 

R Gas  constant 

r Radius 

S Entropy 

T Temperature 

t T i me 

U Rotational  component  of  velocity 

W Velocity  in  the  relative  plane 

z Axial  position 

Ratio  of  specific  heats 
0 Angular  displacement 

\ Second  coefficient  of  viscosity 

H Viscosity 

e Density 

Shearing  stress  vector 
Velocity  potential 
w Rotational  frequency 
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